package rungekutta;

public abstract class RungeKutta {

    protected int n;           // number of the equations in the system
    protected double[] y;      // current values
    protected double x;        // argument
    protected boolean result;

    /**
     * Describes the right part of the differential equation
     * @param x x value
     * @param y y value
     * @param i number of equation in system
     * @return calculated value of the equation
     */
    protected abstract double f(double x, double y, int i);

    /**
     * Processes calculated answer
     * @param x current value of the argument
     * @param y current results
     */
    protected abstract void processAnswer(double x, double[] y);

    /**
     * Solves the system of differential equations
     * @param dx step
     * @param e precision
     */
    protected void solve(double dx, double e, int order, double ae) {
        x = 0;
        double delta;
        double deltaMax = 0;
        double[] ny;
        boolean autodx = false;
        double[] dx1, dx2;
        dx1 = new double[this.n + 1];
        double l, max;
        double linc, ldec;
        if (dx == 0) {
            autodx = true;
            dx = 0.1;
        }
        linc = 0.0;
        ldec = 0.0;
        while (true) {
            deltaMax = 0;
            if (autodx) {
                dx1 = step(dx, x, this.y, order);
                dx2 = step(dx / 2, x, this.y, order);
                dx2 = step(dx / 2, x + dx / 2, dx2, order);
                max = 0;
                for (int i = 1; i <= this.n; i++) {
                    l = Math.abs(dx2[i] - dx1[i]);
                    if (l > max) {
                        max = l;
                    }
                }
                if (max > ae) {
                    if (ldec == dx) {
                        break;
                    }
                    ldec = dx;
                    dx /= 2;
                    continue;
                }
                dx1 = step(dx, x + dx, dx1, order);
                dx2 = step(2 * dx, x, this.y, order);
                max = 0;
                for (int i = 1; i <= this.n; i++) {
                    l = Math.abs(dx2[i] - dx1[i]);
                    if (l > max) {
                        max = l;
                    }
                }
                if (max < ae / 2) {
                    if (linc == dx) {
                        break;
                    }
                    linc = dx;
                    dx *= 2;
                    continue;
                }
            }
            ny = step(dx, x, this.y, order);
            for (int i = 1; i <= this.n; i++) {
                delta = Math.abs(ny[i] - this.y[i]) / dx;
                if (deltaMax < delta) {
                    deltaMax = delta;
                }
                this.y[i] = ny[i];
            }
            if (!result) {
                processAnswer(this.x, this.y);
            }
            x += dx;
            if (deltaMax < e) {
                break;
            }
        }
        if (result) {
            processAnswer(this.x, this.y);
        }
    }

    /**
     * Calculates one step in Runge-Kutta method
     * @param dx step
     * @param x argument
     * @param y current values
     * @order order of the Runge-Kutta method
     * @return new values
     */
    private double[] step(double dx, double x, double[] y, int order) {
        double[] k1 = new double[this.n + 1];
        double[] k2 = new double[this.n + 1];
        double[] k3 = new double[this.n + 1];
        double[] k4 = new double[this.n + 1];
        double[] r = new double[this.n + 1];
        for (int i = 1; i <= this.n; i++) {
            k1[i] = f(x, y[i], i);
        }
        for (int i = 1; i <= this.n; i++) {
            if (order == 4) {
                k2[i] = f(x + dx / 2, y[i] + k1[i] * dx / 2, i);
            } else {
                k2[i] = f(x + dx, y[i] + k1[i] * dx, i);
            }
        }
        if (order == 4) {
            for (int i = 1; i <= this.n; i++) {
                k3[i] = f(x + dx / 2, y[i] + k2[i] * dx / 2, i);
            }
            for (int i = 1; i <= this.n; i++) {
                k4[i] = f(x + dx, y[i] + k3[i] * dx, i);
            }
        }
        for (int i = 1; i <= this.n; i++) {
            if (order == 4) {
                r[i] = y[i] + dx * (k1[i] + 2 * k2[i] + 2 * k3[i] + k4[i]) / 6;
            } else {
                r[i] = y[i] + dx * (k1[i] + k2[i]) / 2;
            }
        }
        return r;
    }

}
